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ABSTRACT 

Radiative transfer and radiation hydrodynamics use the relativistic Boltzmann equation to describe the kinet- 
ics of photons. It is difficult to solve the six-dimensional time-dependent transfer equation unless the problem 
is highly symmetric or in equilibrium. When the radiation field is smooth, it is natural to take angular moments 
of the transfer equation to reduce the degrees of freedom. However, low order moment equations contain terms 
that depend on higher order moments. To close the system of moment equations, approximations are made to 
truncate this hierarchy. Popular closures used in astrophysics include flux limited diffusion and the Mi closure, 
which are rather ad hoc and do not necessarily capture the correct physics. In this paper, we propose a new 
class of closures for radiative transfer and radiation hydrodynamics. We start from a different perspective and 
highlight the consistency of a fully relativistic formalism. We present a generic framework to approximate ra- 
diative transfer based on relativistic Grad's moment method. We then derive a 14-field method that minimizes 
unphysical photon self-interaction. 

Subject headings: radiative transfer — hydrodynamics — relativity 
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1. INTRODUCTION 

Radiative transfer and radiation hydrodynamics use the 
ultra-relativistic Boltzmann transport equation to describe the 
kinetics of photon. The radiative intensity, which is pro- 
portional to the photon distribution function, is a seven- 
dimensional hypersurface embe dded in eight-dimensional 
phase space (Misne r et al.l 11973b . Unless the radiative in- 
tensity is in equilibrium, or the problem is highly sym- 
metric, it is difficult to solve the radiative transfer equa- 
tion either a nalytically or numerically (see standard text- 
book such as|^ andrasek IiarJ|1960l;lMihalas & Mihalaslll984t 
Rvbicki & Lightman Tl986t lShulll991t lPeraiahll200U iCastori 
2004; Graziani 2006). On the one hand, there is no successful 



theory to reduce the complexity of radiative transfer. On the 
other hand, although numerical algorithms such as direct dis- 
cretization of the radiative transfer equation and Monte Carlo 
methods exist, they are computationally too expensive. 

In astrophysics, very often we are interested in the radiative 
energy and flux instead of the intensity. If the radiation field is 
smooth (in terms of directions), we simply take angular mo- 
ments of the radiative transfer equation and solve only for the 
frequency-dependent moment equations. The zeroth-, first-, 
and second order angular moments of the intensity carry clear 
physical meanings. They are the radiative energy density, ra- 
diative flux, and radiative stress tensor, respectively. For many 
classic problems in astrophysics like stellar atmospheres, the 
global symmetry of the system is used to further reduce the 
degrees of freedom. Analytical or numerical solutions are 
then obtained by solving the reduced frequency-dependent 
moment equations. This approach is naturally extended to 
radiation hydrodynamics, which not only describes how the 
(movin g) media radiates, but also how radiation feeds back to 
matter dMihalas & Mihalasll984tlM"ihalasl200UlCastorj2004t 
Kj-uiTiholz_et_alJ|2007t). 

The dynamic equations of the radiative stress tensor con- 
tain terms that are related to the third order moment of 
the intensity, while the equations of the third order mo- 
ments depend on the fourth order moments and so on. 

ckchan@cfa.harvard.edu 



In order to close the system, we need to make approx- 
imations and truncate the moment hierarchy. This is 
known as the closure problem in radiative transfer. Pop- 
ular closure schemes in astrophysi cs include variant forms 
of flux limited diffusion (F LD. | Chang& Cooper! 119701: 

01 Pomranind l 19831; iLevermore 
P/v ap proximations (some ti mes with diffusion cor- 



rections, see lOlson et ail EOOOl ISeibold & Frank! 120091: 
iMcClarren et all l2008t iRavishankar et all 120101 "and refer- 
ence therein ), the Mi closure (equivalent to max i mal entropy 
closu re, see Uanka et all 119921 iSmit et all I2000t iFarris et all 
20081). and variab le Eddington factors (VEF.lFeautrierlll964 



Pomraningl 119691: |Auer & MihalaJ 119701: iFukud 12008a 1 



2009fl They have been successfully applied to many astro- 
physics problems, yet their underlying assumptions (what is 
the optimal form of a flux limiter for a particular problem? 
why is the entropy of photons maximized? how many mo- 
ments do we need?) are rather ad hoc and do not guarantee to 
capture the correct physics. 

In hydrodynamics, the Chapman-Enskog theory expands 
the particle distribution in power series of the Knudsen num- 
ber, which is the ratio between the particle mean free path 
and the typical length scale in the problem. The zeroth 
order expansion gives rise to ideal hydrodynamics, while 
the first order terms allow us to calculate transport coeffi- 
cien ts such as kinematic viscosity and thermal conductiv- 
ity dHuandll965tlC hapman & Cowling! [l970t iLibol fl979l 
Cerci gnani & K remer 20Q2). We would like to apply a simi- 
lar technique to radiative transfer in order to solve the closure 
problem. Unfortunately, photons do not self-interact. The 
photon mean free path is determined by the extinction coef- 
ficient, which is a material property. In the cases when the 
photon density is high but the media is optically thin, radia- 
tion dominates the energy budget but cannot be thermalized. 
The Chapman-Enskog series for photon distribution function 

1 The term VEF has recently been abused in the literature. The original 
Auer & Mihalas ( 1970) algorithm uses an iterative solver for the Eddington 
factors, which is technically not a closure relation. The closures such as M i 
should really be called moment methods with prescribed Eddington factors 
(Fu 1987; Park 1990; Koerner & Janka 1992). 
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fails to converse (Mihalas 2001). 

In addition to the closure problem, as studies of astro- 
physical objects become more detailed, complicated sub- 
structures such as convection in stars and turbulence in ac- 
cretion disks are always found. The radiating medium moves 
non-uniformly, which makes the computation of spectra more 
complicated due to the spatial-dependent Doppler effect. The 
standard approach is to Taylor expand the extinction and 
emission coefficients and to derive moment equations to at 
least O(v/o), where v and c are the speed of the media 
and the speed of light. The resulting equations have terms 
that are physically impor tant even in the limit v/c — > 
(Mihalas & Mihalas 1984). Keeping track of these higher or- 
der terms in different p hysical regimes is a challenging task 
dKrumholz et al.l [20071) . Moreover, there are astrophysical 
systems such as the inner regions of accretion disks and ultra- 
relativistic jets that move with relativistic speed so the de- 
scribed approach converges too slowly. 

If radiative feedback is weak, we can solve the hydrody- 
namic equations and compute the radiative spectra us ing post- 
processing (see iNoble et all 120071; IChan et alJ 120091) . How- 
ever, there is a large set of problems for which the feedback is 
important. We have to solve the frequency-dependent three- 
dimensional moment equations at every single time step. The 
full system of radiation hydrodynamics is highly non-linear. 
Studies of turbulence show that reducing the number of spatial 
dimensions and enforcing symmetr y in these systems can p ro- 
duce fundamentally wrong results ( Kraichn anl 1 9671 [1971). 

To summarize, we have listed three major difficulties in 
using moment methods for radiative transfer and radiation hy- 
drodynamics: 

(i) Because of the collisionless nature of photons, standard 
methods in kinetic theory such as the Chapman-Enskog 
theory fails to provide a physical closure. 

(ii) When the non-uniform motion of media is considered, 
it is difficult to take properly the spatial-dependent 
Doppler effect into account. 

(iii) If the radiative feedback is important, the naive ap- 
proach is computationally too expensive to couple ra- 
diation back to hydrodynamics. 

Hydrodynamics depends only on the frequency- 
integrated radiative force. Many existing codes, 
therefore, reduce the degrees of freedom and resolve 
problem (i i i) by using frequency-integrated equations 
dStone et al.l 119921 iTurner & Stond I200U [Haves & Normanl 



20031 iGonzalez et alJ 120071 IKrumholz et alJ 120071 
Aubert & Teyssieri 120081: iFarris et all 120081) . If the radiative 



spectrum is needed, it is always possible to post-process the 
numerical solution. 

Once we integrate over frequency, the moment equations 
are fully relativistic. The radiative energy density, momen- 
tum, and stress tensors together form a covariant stress- 
energy tensor. We can then ignore the complication of 
mixed frame formalism. This relativistic approach is origi- 
nall y developed to solve radiative transfer in genera l relativ- 
ity (Anderson & Spiegel 1972; Schmid-Burgk 1978; Thorne 

2 Another fundamental difficulty of radiation hydrodynamics comes from 
the large separation between the radiative time scale and the (hydro-) dynamic 
time scale. Nevertheless, this difficulty raises not because of the moment 
method. It is a generic property of non-relativistic astrophysical systems. We 
will therefore leave it out from this paper. 



[19811 lUdev& Israeli fl982l) . It is well adopted in numeri- 
cal general relativistic hydrodynamics dAguirre et al.l 120051 
lParkll2006h lTakahashi 2007: IFarris et~aT] l2008). However, we 
believe that its true advantage is in the covariant equations, 
w hich t r ivially addresses problem (ii). 

iGradl (| 19491) proposed a moment method which expands the 
ratio between the particle distribution and the equilibrium dis- 
tribution in a multi-dimensional polynomials of the momen- 
tum. Instead, it treats the moments as fundamental fields and 
keeps track of their evolutions. This method is independent of 
the mean free path. Once the reference frame and temperature 
are chosen, Grad's coefficients are simply linear transforms of 
the momentum moments, which converge exponentially fast 
when the distribution f unctio n is well behavecQ. Because pho- 
ton transport is linear, IGradT s moments method provides the 
most natural way to res olve pr oblem (i). 

The idea of applying Grad's moment method to stu dy ra- 
diation is not new. A classic paper by iThornd (119811) pre- 
sented very detailed moment formalisms for relativistic radia- 
tive transfer by using projected, sy mmetric, trace-free ten sors 
dThorndll980l) . Shortly after that, [Udev & Israel d!982l) de- 
rived a 14-field approximation for general relativistic radia- 
tive transfer based on IGradT s method. However, as far as we 
know, thes e attempts were only used in one-dimensional prob- 
lems (e.g. Zam pieri et al.lfT996h . On the one hand, this is pos- 
sibly due to the general concept that Newtonian formalisms 
are always simpler than relativisti c formalisms. On the other 
hand, IThornd and lUdev & Israeli s 14-field methods are in- 
de ed complicated compared t o the standard 0(v/c) equations 
in iMih alas & Mihalai dl984l) . Their advantages only appear 
when we consider more physical regimes such as the ones de- 
scribed in IKrumholz et al.l (12007b . 

Based on the above points and some additional physical un- 
derstanding of moment methods, we propose a new class of 
closure schemes for radiative transfer and radiation hydrody- 
namics in this paper. The paper is organized as the following. 
In the next section, we introduce our notations and review the 
standard radiation hydrodynamics equations. In section[3] we 
describe [Gradf s moment method and its ultra-relativistic gen- 
eralization. We also summarize the linear transformations that 
we derive in the appendixes. In section @] we study carefully 
ideal hydrodynamics and some existing closures for radiative 
transfer. They provide us important physical insight, which 
leads to a new class of closure schemes. Although the scheme 
is generic, we specifically look at the 14-field approximations 
and derive the closure equations in section[5] Finally, we con- 
clude this paper in section [6] 

2. NOTATIONS AND STANDARD EQUATIONS 

We use the component notation in Mis ner et al.1 (119731) 
though out this paper. Greek indices run from to 3 and Latin 
indices run only from 1 to 3. The metric tensor is denoted 
by g a ° and the metric signature is (—,+,+,+). A point in 
spacetime is denoted by x a = (ct, x l ) with c being the speed 
of light. The Einstein summation convention is used unless 
specified otherwise. 

We start from the standard radiative transfer equation. Let v 
be the photon frequency, n % be a three-dimensional unit vec- 
tor. We define l v = I(t,x l ;i , ,n t ) be the specific radiative 
intensity, r\ v be the total emission coefficient, and Xv be the 
total extinction coefficient. The intensity is then governed by 

3 It is a simple application of Darboux's principle. See, for example, Boyd 
1 1999), on convergent properties of series expansions. 
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the radiative transfer equation 



1 d i d w 



(1) 



We first describe the standard moment method s 
dChandrasekhari [l960t iMihalasI fl97l IMihalas & Mihala 
1984 lRvbickT & Lightmanl 1 19861: IShul [T99 It iPeraiahl |200T 
Castor! 120041) . Taking the zeroth and first angular mo 
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mentum of the whole transfer equation (fTJ, we obtain two 
frequency-dependent moment equations, 



d t E v + a i Ft-- 
d t Fi + c 2 d,py : 



-cG° 
-c 2 Gt 



(2) 
(3) 



The quantities on the left are the frequency-dependent radia- 
tive energy density, radiative flux, and radiative stress tensor, 



dfll v , 



F l v = JdV.I„n l , 



P*> = I (IV l.n' ir'. 



(4) 



(5) 



(6) 



The ones on the right are the radiative energy and moment 
inputs (to the radiating medium), 



(7) 
(8) 



Note that the four-tuple (E v ,Fl") is not a four- vector so 
equations (O and (01 do not form a covariant equation. The 
quantities E v , F % u , and P % J , as well as G® and G l v , all depend 
on the reference frame. This is exactly the difficulty in solving 
angular moment equations when the radiating medium moves 
non-uniformly — we cannot Lorentz transform (E V ,F*) to 
obtain their values in a different frame. 

We follow IMihalas & Mihala! (11984 to derive the covari- 
ant formulas as we suggested in the introduction. Let h be 
Planck's constant and n a = (1, ri 1 ) be the "unit" null vector. 
The photon four-momentum is p a = {hv/c)n a . The Lorentz 
invariant photon distribution function / is related to the inten- 
sity by the equation 



J(i,x; (9) 



We can now rewrite equation (03 as a ultra-relativistic Boltz- 
mann transport equation, 



P7 ;Q = e-xf, 



(10) 



where the subscript - a denotes covariant derivative. The 
Lorentz invariant emission coefficient 



e(i/) = 



and extinction coefficient 



hv 

x{v) = — Xv 
c 

are usually assumed independent of /. 



(11) 



(12) 



Using the covariant volume element, d 3 p/p°, we can take 
four-momentum moments of equation ( fTOb to arbitrary order, 

(13) 



ai + l 



-G° 



The moment tensor on the left hand side is defined by 

Raia2 -a l + 1 _ c f fPf ?aipa2 . . o 

J p" 

while on the right hand side, we have 

' d 3 p 



G° 



(14) 



(15) 



We refer G axa% "' ai as the moment extinction term because 
— G" 1 " 2 " is usually calle d the moment production term 
(ICercignani & Kremerll200"2h . 

It is easy to connect the relativistic formulas with the non- 
relativistic ones. Using the first order equation, the radiative 
stress-energy four-tensor is 



= - I dl'dV. /..;,••// 



1 



E Fi/c 
F l /c P l] 



(16) 



The frame-dependent components are the frequency- 
integrated radiative energy density, radiative flux, and 
radiative stress tensor. Similarly, the radiative four-force is 



G a = X - \(hnlV.[\ i :i„ - ,/„),/" = 



G° 
G" 



(17) 



The temporal and spatial components are the frequency- 
integrated radiative energy and momentum inputs. We remark 
that higher angular moments of the intensity, even though they 
are frequency-integrated, do not have counterparts in equa- 
tion (fT3ll. 

For completeness, we also write down the equations for rel- 
ativistic hydrodynamics. We use v % to denote the fluid veloc- 
ity. The four- velocity is u a — j(c,v l ), where 7 is the Lorentz 
factor. The material four-momentum is given by 



rpa 



pu" 



(18) 



where p is material density in the local Lorentz rest frame. 
For simplicity, we assume perfect fluid so that the material 
stress-energy tensor takes the simple form 



T af} = P hu a uP +pg afi . 



(19) 



Here, h = e/c 2 + p/pc 2 is the specific enthalpy, e is the 
specific internal energy density including the rest energy, and 
p is the thermal pressure. 

After all the definitions, standard radiation hydrodynamics 
can be summarized in three tensor equations, namely, the con- 
tinuity equation, 

T% = 0, (20) 



the covariant Euler equation, 

T aP ,p = G a , 
and the radiative stress-energy equation 

R afi ,p = -G a . 



(21) 



(22) 



The above equations are Lorentz covariant. Nevertheless, 
we refer them as "lab frame equations" in order to distin- 
guish from the radiation fiducial frame and the fluid comoving 
frame. 
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3. GRAD'S MOMENT METHOD 

Because of the perfect fluid assumption, the material part of 
standard radiation hydrodynamics is closed by equation dl9l l 
and an equation of state. We will refer this as the ideal (rel- 
ativistic) hydrodynamic closure. Conversely, there are nine 
unknowns in the trace-free symmetric radiative stress-energy 
tensor but only four independent equations. The transport part 
of radiation is subject to the closure problem. In addition, 
it is unclear how to rel ate the radi ative four-force to other 
macroscopic quantities. Grad (1949) proposed evolving the 
moments as fundamental fields and use them to reconstruct 
the particle distribution function. The unknown components 
of the moments and the moment extinction terms can then be 
evalu ated self-consistently. In this section, we will generalize 
Grad's moment method to work with radiation. 



3.1. Grad's Expansion 

We generalize Grad's moment method in a covariant form 
and expand / as the follow multi-dimensional power series 

/ = / (0) 



a a p 



at each point x a . In the above expansion, we choose 

2/h 3 



f 



(o) = 



exp(-p a U a /8) - 1 



(23) 



(24) 



be the equilibrium (black body) photon distribution function, 
where U a is a four-velocity of some fiducial observer and 8 
is the energy scale for the equilibrium distribution. The co- 
efficients a, a a , a a /3, ... are tot ally symmetric with trace-free 
spatial parts (see appendix[A]or Thorne 1983). 

In the fiducial reference frame, ~p a U a reduces to the pho- 
ton energy hv. Using the energy scale 8, we define a dimen- 
sionless quantity 

£ = hv/0. (25) 

The equilibrium distribution reduces to f(°> = (2/h 3 )w(^), 
where 

w & = J7T^ (26) 

can b e interpreted as a dimensionless weight. We can rewrite 
iGradT s expansion as 



/=£«,(0(a 



+ £a a n a + £ 2 a a pn a nP + 



(27) 



The rescaled coefficients a, a a , a Q ^, ... are all dimensionless. 

Note that the weight iu(£) depends explicitly on direction 
after a Lorentz transformation — this is simply the Doppler 
effect. The coefficients in expansions d23l and d27l ) are spe- 
cific for the (at the moment, arbitrarily) chosen energy scale 
and fiducial frame. Nevertheless, the infinite series contains 
full information of the original Lorentz invariant distribution 
/. We can Taylor expand t he an isotropic part of the weight 
and obtain another set of iGradf s coefficients in the other 
frame. The distribution is fiducial frame dependent only after 
we truncate the expansion at certain order. 

3.2. Solving for Grad's Coefficients 

In principle, we could deri ve evo lution equations for lGradf s 
coefficients by substituting IGradf s expansion (|27T i into the 
ultra-relativistic Boltzmann equation (ITOb but w e wou l d then 
lose the physical intuition. Instead, we follow iGralCl 949) 
and treat the moments as fundamental fields. The left hand 



sides of the moment equations ( fT3T > suggests that we need to 
keep track of R a i°'2— a i (components that have at least one 
in the indices). The flux terms R ni2 (components with 
no in the indices) and the extinction terms Q a i a 2 - °ii m ust 
be solved in terms of them. 

Taking moment is a linear operation. We can skip the distri- 
bution function / and directly write down a linear transforma- 
tion between the coefficients and the moments. Let the Euler 
script 

,3 / -x '-i 



X J l c 



R 



OL\OL2 - - - CK( 



(28) 



be the dimensionless moments, where A = he/ 8 is the photon 
mean separation. Using appendix iBl the relation between the 
dynamic variables and Grad's coefficients can be summarized 
in the following hierarchy of matrix equations. 



(29) 



R° - 




r W 2 


w 3 


W 4 




' a 


£00 




w 3 


W 4 


w 5 ■•• 




a 


£000 


= 8tt 




w 5 


We ■■■ 




a 00 



ftOi 



W 4 2W 5 
W 5 2W 6 



a, 

«0i 





1 




~3 



^000 



16tt 
~15 



W 6 



(30) 



(31) 



where Wi is the shorthand of the integral Wi = J °° dt;w({;)£, 1 . 
Substituting our weighting function (|26| |. the integral has the 
closed from solution Wi = T (l + 1)C0 + !)• Some numerical 
values are listed in equation (IB4t . 

Of course, there are infinitely many equations in the hierar- 
chy. We just present the ones that are useful for this paper. It 
is easy to solve the coefficients by inverting the above trans- 
forms. We will provide some examples in section |4] and 

3.3. Computing the Flux Terms 

Once we obtain IGradf s coe fficients, we can use them to 
compute the flux terms. Grad's moment method is linear in 
the fiducial frame. The linearity naturally form a class of clo- 
sure schemes. Since we fix the weight the only free- 
doms in the closures are the energy scale 8 and the fiducial 
reference frame corresponds to U a . 

We will discuss how to choose the fiducial velocity U a in 
section[4] For now, we use the subscript r to indicate the radi- 
ation fiducial frame, 



■■■A?' a R 



(32) 



where A? a is the Lorentz transformation matrix. We use the 



l r/3 



result derived in appendixlBlto compute the flux terms. 

&7T / 



TJ = 



3 

16-7T / 

~15 V 



W?,a 



3W f 



35 V ryfc 



6 a r0ij 

5 



(33) 



(34) 



(35)' 
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To obtain the flux terms in the lab frame, we simply apply an 
inverse Lorentz transformation to R" lCt2 '" ai . 

Note that the dimensionless photon number density 9?° does 
not enter the closure equations, we can use it to choose the 
energy scale. We require 9 to approach fceT in local thermal 
dynamic equilibrium, i.e., 3?° = 8ttW2, which leads to 



h 3 ( 



8nW 2 



1/3 



(36) 



3.4. Computing the Extinction Terms 

Photon emission and absorption are fluid properties. They 
take their simplest form in the fluid comoving frame. There- 
fore, we use the fluid temperature k^T and the fluid four- 
velocity u a to replace 9 and U a in Grad's expansion (see ap- 
pendix 0. 

To compute the extinction terms, we apply Lorentz trans- 
formation to boost the moments to the fluid comoving frame, 



a i a. 2 ■ ••ot-i 



A? 1 A Q2 



(37) 



We follow the same procedure described in section I3T2I to ob- 
tain the comoving coefficients af, a,f a , a{ a p, ... The subscript 
f here indicates the fluid comoving frame. Introducing the 
notation 

6 w 



a. i a. 2 • "C^z 



and letting the Euler script 

\4 / \ J 

A I C 



\k B T 



Q a l a '2 '"Gil 



(38) 



(39) 



be the dimensionless comoving extinction terms, the deriva- 
tion in appendix|C]gives us the following linear transforms. 



St 

s?° 





\x 2 


x 3 






x 4 


8tt 


x 4 


x 5 



5? 



X4 2X5 
X5 2Xq 



b f 
b{ 
h 00 



bit 
bfoi 



(40) 



" 5? ' 


1 




ij 167T 


'x 6 




bfij 




~3 




T HT 









(41) 



(42) 



There is no freedom left in the above equations because T and 
u a are both fixed as fluid properties. We can apply inverse 
Lorentz transform to boost the extinction terms back to the 
lab frame. 

4. PHYSICAL IMPLICATIONS IN CLOSURE SCHEMES 

It seems that we have all the equations to close the radiative 
moment hierarchy. Unfortunately, iGradT s moment method 
has some arbitrariness in the choice s of t he fiducial frame ve- 
locity U a as we remarked in section [3731 In order to get some 
physical insights to constrain these arbitrariness, we need to 
understand closure schemes at a more fundamental level. 



4.1. The Origin of Non-linearity 

Let us think deeper about the very successful ideal hydro- 
dynamic closure. It only requites the left hand side of the 
Boltzmann transport equation, which is linear, to derive the 
Euler equation. Taking moments is also a linear operations. 
So why is the resulting equation non-linear? 

Without loss of generality, we use non-relativistic hydrody- 
namics for this discussion. Let m be the particle mass, 



- / -/■'>./>' 
m 



(43) 



be the fluid velocity, and p n = p l — mv l be the particle mo- 
mentum in the comoving frame. The particle distribution / 
is always well approximated by a Maxwellian in the fluid co- 
moving frame because of collisions. We can choose a linear 
closure in the comoving frame, such as p — (2/3) pe, yet the 
non-linear inertial term will always appear in the lab frame, 

— / d 3 pfp l p 3 =muV [ d 3 pf + — { d 3 pfp H p' 3 
m J J m J 

= pv i v j +pg lj ■ (44) 

Therefore, non-linearity does not just come from the clo- 
sure approximation. It is a direct consequence of Galilean 
transforming the second (or higher order) moment with a ve- 
locity that depends on the distribution. In other words, hy- 
drodynamics is non-linear because of Galilean symmetry and 
fluid particle self-interaction. 

For radiation, we replace Galilean symmetry by Lorentz 
symmetry in the above reasoning. Photons do not self-interact 
so they can only be thermalized by external medium. In the 
diffusion regime, a few moments can fully describe the pho- 
ton distribution function. It does not matter rather we use the 
fluid velocity u a or some radiative quantities to solve for the 
radiation fiducial velocity U a . Both choices can lead to fully 
relativistic and well behave closures. 

In the free-streaming regime, however, the radiating 
medium does not contribute much. On the one hand, setting 
jja _ u a j s non . sense On the other hand, employing the lab 
frame for linearity causes closures, such as the P/v closures, 
violates the Lorentz symmetry. We believe that Lorentz sym- 
metry is an important property so we ensure that it is satis- 
fied in the proposed scheme. 

4.2. Unphysical Photon Self- Interaction 

It is educational to study the simplest Lorentz invariant clo- 
sure. We need at least one variable for the energy scale and 
three variables for the fiducial velocity. By symmetry, the ra- 
diative stress-energy tensor must take the form 



R? = 



3P 














p 














p 














p 



(45) 



4 There is an exception. If we need to implement a radiative transfer solver 
on a low resolution grid, the spatial discretization introduces fiducial frame 
dependence because of truncation error. It this case, it may not be a bad 
idea to just give up Lorentz invariance and use high order linear closures. In 
some sense, this is the philosophy behind Latt i ce Boltzmann methods (see, 
for example. ISucci et aU1 993; He & Luo 1997; Lallemand & Luol200q) . 
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in the fiducial frame. Boosting it back to the lab frame, we 
obtain the standard stress-energy tensor for photon fluid, 



R a0 = PJJCXJJ0 +P g^ 



(46) 



Recalling equation (fT6*b and compare different components, 
we obtain 



e | — u°u° - 1 1 r>. 



F l = - U U l P, 
c 

/"- (4'"' ' • ) < p 



(47) 
(48) 
(49) 



Let us treat the radiative energy E and flux F l as the funda- 
mental fields. The closure problem reduces to solving P %3 in 
terms of E and F l . The solution is simply 



py = h Pq 13 

c 2 {E + P) y 



with 



E / / 3F l Fi 



(50) 



(51) 



The above equations describe exactly the Mi closure al- 
though the deriva t ion is different f rom the standard ones 
(IJanka et al.1 [l99l ISmit et all l2000t iFarris et alJ l2008h . In 
fact, our derivation is probably more general since the only 
assumes are the number of fundamental fields, isotropy in the 
fiducial frame, and Lorentz invariance in the closure. There is 
no assumption about the spectrum. 

The above derivation is fully relativistic so the information 
propagation speed is limited by the speed of light. The closure 
appears to have correct diffusion and free-streaming limits. 
However, the derivation implies that the photons always ther- 
malize themselves even without interacting with the radiating 
mediurrQ. We can interpret the truncation errors as unphysical 
photon self-interaction. Our closure problem now reduces to 
minimizing this unphysical effect. 

4.3. Moment Decomposition 

There is no freedom left in the Mi closure. We have to 
introduce more dynamic variables to reduce the unphysical 
photon self-interaction. The simplest trial is to include the 
zeroth moment of the radiative transfer equation, 



R - a — G. 



(52) 



The above equation and equation d22l form a 5-field method. 
Compare to M\, the extra information governed by equa- 
tion d52l allows us to go beyond the gray approximat ion. It 
is easy to derive the extinction terms by using section [341 or 
appendix [C] For the purpose of this paper, we only focus at 
the closure problem. 

The structure of the 5-field method is similar to relativis- 
tic hydrodynamics (|20| i and (BTT i. where the fluid velocity is 
equal to the fiducial velocity. I.e., the fiducial frame is chosen 
so that the particle density flux vanishes. The heat flux, vis- 
cosity, etc are then defined in such a frame. There is a special 

5 Recalling that the M\ closure is originally developed to maximize the 
entropy (Janka et al. 1992). 



name associated with thi s choice: the | Eckartl (1 19401) decom- 
position. Alternatively, lLandau & Lifshitzl dl959l) proposed 
shifting the fiducial velocity so that the heat flux vanishes in 
the fiducial fram e, = 0, where the subscript l denotes the 

lLanda u-Lifshitzl frame. 

The fiducial velocity in the Eckart decomposition describes 
the particle flow, while in the lLanda u-Lifshitz decomposition 
it describes the energy flow. When the radiation is not in local 
thermodynamic equilibrium, the two decompositions result 
different fiducial frames and have different level of unphys- 
ical photon self-interactions. However, the 5-field method is 
too restrict to describe non-equilibrium effects. The photon 
density flux R 1 is parallel to the photon energy flux R° l in the 
fiducial frame, 

K = *?. (=3) 

The two decompositions are therefore identical. We must use 
higher order moment methods to reduce unphysical photon 
self-interaction. 

5. A PHYSICALLY MOTIVATED 14-FIELD METHOD 

To describe non-equilibrium effects, we follow|Gri|dl949) 
to truncate the expansion (|27T i at se cond o rder. The photons 
are described by an ultra-relativistic lGrad 's distribution func- 
tion, 



f 



(2) 



h 3 



^a a n c 



/• (54) 



Because a a /3 is symmetric and has trace-free spatial parts, 
there are only nine independent components. Taking a and 
a a into account, the polynomial has fourteen independent co- 
efficients. They can be solved by the fourteen time dependent 
fundamental fields R°, R 0a , and R 0al3 . 

To derive our 14-field method, we first reduce the flux equa- 
tions to include only the non-vanishing coefficients, 



2W 4 a r 



X 3 = 



T I 



(55) 
(56) 
(57) 



The dimensionless moments in the above equations are 
rescaled by the energy scale 8 chosen in equation ( T3"6*i >. Note 
that a, ao, and aoo do not appear in the flux equations. There- 
fore, we can skip equation d29l and just inverse the reduced 
forms of equations ( f30b and QTl ), 



(58) 



mOi 
Jv r 


Stt 


' 14^4 


2W 5 ~ 




al 






w 5 


2W 6 




a? 



Eliminating lGradT s coefficients, we have 



mi 
J\.. 



W 3 W 6 - W 4 W 5 
WiWa - W 5 W 5 



rpU 



w 3 w 5 - W 4 W 4 
W 4 W 6 - W 5 W 5 ~ 



mo 



W e V 3 1 T 



(59) 



(60) 



(61) 



Equation (l60l l. ( f6Tb . and d57l i almost form a closure relation. If 
we fixed the fiducial frame to the lab frame, it would become 
a linear (but frame dependent) closure similar to P%. 
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Equation ( f60b shows that the lEckartl and lLand au-Lifshitz 
decompositions result different fiducial frames. In addition, 
we can introduce a third order decomposition so that the fidu- 
cial frame is chosen to satisfies the requirement 3J°° 4 = 0. 
Equation d60t becomes 



X = 
X. = 



0.10353?°°* in me |Eckarj frame, (62) 
-0.05483?°°' in the lLandau-Lifshitzl frame, and (63) 



0.52993?" 1 in the third order frame, 



(64) 



respectively. We can refer the lEckartl and Landau-L ifshitzl 
frames as first and second order frames. 

The sign difference between 3?* and 3?° 0i in eq uation d63l 
has interesting meaning. It indicates that the Landau-Lifshitz 
fiducial velocity lies between the other two. Appendix [D] 
demonstrates that higher order velocities are less sensitive to 
the distribution function. The non-linear terms have small 
contribution to the overall dynamics. Therefore, we conjec- 
ture that higher order decompositions introduce less unphysi- 
cal photon-self interaction. In order words, we propose using 
the third order decomposition for the 14-field method. 

6. CONCLUSIONS 

In this paper, we propose a class of physically motivated 
closures for radiative transfer and radiation hydrodynamics. 
We start by showing the advantages of frequency-integrated 
schemes, which are fully relativistic. The transport terms and 
extinction terms can be easily evaluated in different reference 
frames. We then apply Grad's moment method to compute 
the flux terms and the extinction terms from the fundamental 
fields. The truncated Grad's series is energy scale and fiducial 
frame dependent. For the extinction terms, the fluid comoving 
frame and temperature are the natural choices. For the flux 
terms, however, we propose using high order decomposition 



to reduce unphysical photon self-interaction as well as using 
photon number density to obtain the energy scale. 

We believe that this paper clarifies some implicit assump- 
tions in standard closures and points out the importance of 
fiducial frame. Although we only present the 14-field method, 
it is straightforward to derive arbitrarily high order closures 
from our formulas (see the Appendixes). With minimal modi- 
fications, our closures are also applicable to neutrino transport 
and relativistic rarefied gas dynamics. We expect the methods 
derived from our new framework to outperform existing ones. 

Of course, there are many open questions associated with 
the proposed schemes such as limiting behavior and linear 
stability of the theory. We will address these issues in subse- 
quent papers. We are also implementing the 14-field method 
and integrating it with some hydrodynamic solvers. We will 
perform verifications and validations before applying the al- 
gorithms to study astrophysical systems. As a final remark, 
moment method s are not optimal for solving all problems 
dGrazianil [2006). For example, for situations with strong 
beaming, we are be tter off if we use other techniques such 
as ray tracing (see 



Finlator et al. 2009 



Abel & Wa ndelt 200^ iTrac & Ced 120071 



and reference therein). 
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APPENDIX 

A. GRAD'S EXPANSION AND SYMMETRIC TRACE-FREE TENSORS 

Following the notations in section [2] we use p a = p°(l, n l ) to denote the four- momentum of a massless particle, where n l is 
some spatial unit three-vector. We assume there exists energy scale 9 to describe the particles in some fiducial reference frame. 
Note that 9 is not necessary related to the temperature. Let U a be the four-velocity of the fiducial frame, we can define £ as a 
dimensionless measurement of energy for each particle, £ = —p a U a /0. For photon with frequency v, it reduces to £ = hv/0'm 
the fiducial frame, where h is Planck's constant. 

Considering a fixed position x a in the fiducial frame, we can expand the local particle distribution function in a weighted 
multi-dimensional polynomial of p a , 



+ oft/ 1 + aftft/ 1 / 2 + fift^/W 3 



(Al) 



In the above equation, g s is the number of available states, which equals 2 for "photon gas". The function w(£) is a dimensionless 
weight. When it is well chosen, the integral for each term in the polynomial is guaranteed to converge when we take moments. 
Bear in mind that £ is subject to Doppler shift. Therefore, the above expansion depends on the choice of the fiducial frame, i.e., 
U a , in addition to 9. The coefficient <if) 1 B 3 ...B l is a i-th rank tensor, which has dimension (c/6) , It is implicit in the notations 
that these coeffic ients a , Oft, Oftft, ... depend on the position x a . 

In the original iGradP s moment method for non-relativistic and non-degenerated gas, because of the Maxwellian distribution 
in equilibrium, the weighting function to(£) is chosen to be a multi-dimensional Gaussian, which turns out to be the weight 
for the orthogonal conditions of Hermite polynomials, resulting a Hermite series expansion. By the same token, a natural 
weight for relativistic non-degenerated gas is an exponentially decaying function because of the Maxwell- Jiittner distribution. It 
corresponds to the orthogonal condition of Laguerre polynomials. Using this choice, 9/k-Q converges to temperature in the limit 
of local thermodynamic equilibrium. In this paper, however, we use w(£) = [exp(£) — 1] _1 for radiation so the polynomial 
expansion reduces to unity in the limit of black body radiatiorfl It is not a weight for classical/standard orthogonal polynomials. 
Nevertheless, because of its asymptotic behaviors, it is easy to verify that all the useful integrals are guaranteed to converge. 



6 The moment method proposed here is very general. When applying to Wl = (1 - 2~»)r(l + lYU + 1) instead of equation fig), 

neutrino transport, we can chose w{£) = [exp(§) + 1] . This results 
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From expansion ( lAU . it is clear that the tensor ap x 2 ...^ is totally symmetry. In addition, it has trace-free spatial parts because 
p a has only three independent components. One way to verify this is to rewrite the above expansion as 

/= p w(0 (a + £a Pl n^ + ^a^n^ + e^ M n^n 02 n 03 + • • •) (A2) 
= || w(0 (a + £a + £ 2 a 00 + fa 00 o + ■••)+ ^(a h + 2£a 0il + 3£ 2 a 00 ii + • • -)n ,;i + 



(A3) 



where the rescaled coefficients a, ap ± , a^ 1( g 2 , ... are all dimensionless. We can now employ the results in iThornel {1980) to 
conclude a^...^ must be a linear combination o f sy mmetric trace-free tensors. Therefore, ap 1 p a ...p l has trace-free spatial parts. 
It is easy to see the close relationship between (lAlb and hydrogen energy eigenstates by using the identity between spherical 
harmonics and the basis set of symmetric trace-free tensors. 

Counting degrees of freedom is another way to verify the trace-free property. Recalling the zero moment of the ultra-relativistic 
Boltzmann equation yields one dynamic equation; the first moments yield four dynamic equations. However, the second moments 
yield only nine equations because photons are massless. In order to match the degrees of freedom at each order, a / 3 1 p 2 ...p l has to 
be totally symmetric with trace-free spatial parts. 

B. MOMENTUM MOMENTS OF MASSLESS PARTICLES 

We are now ready to derive the equations for computing four-momentum moments, iJ" 1 " 2 '"" 1 , from an ultra-relativistic one- 
particle distribution function /. First, we separate the ^-dependence and angular-dependence of the particle moment p a by 
writing 

d 3 v 

— K-fP P ■ ■ -P 
p" 

p"dp" dQfp ai p a2 ---p°" 

1+2 

C 



d£dtlft l+1 n ai n a2 



1+2 

She'- 
ll 3 



- d^dnw^)(^a + ^a^n Pl +^ 2 a M2 n^n^ 2 + ---y i+1 n ai n a2 ---n a '. (Bl) 

We are free to rearrange the indices because iJ a i°a - " a i is totally symmetric. Let m be the total number of non-zero indices, 
without loss of generality, the moment can be written in the form Jl°—°m2—i m p ac k m g a u me zeros to the beginning and 
leave the spatial indices at the end. We can separate all the ^-dependence and angular dependence, 



Totally I indices 

j^aia 2 ---ai j^O • • • O11I2 • • • im 

1+2 



(B2) 



= |f ( c) / ^ dQ Wl+1 (" + ^ af3inPl + t 2 ^!^ 1 ^ 2 + ■■■) n %1 n t2 ■ ■ ■ n % ™ 
-J d(dQ [wi +1 (a + £a + e 2 a o + ■ ■ •) n^n 12 ■ ■ ■ + 
wi+2 (a h + 2 £a Qjl + 3 fa o n +•••") n n n 12 ■ ■ ■ n lm n 31 + 

Wi +3 (a jlj2 +SS,a Qnj2 +6^ 2 a o 3l i 2 H ^n ll n 12 ■■■n %m n n n n + 

wi+i(a n j 2 j :i +4^a j u - 2 j 3 + 10 ,,, H ^jn ll n 12 ■ ■ ■ri tm n n n 32 n : ' 3 + ■ 

where we have introduced the shorthand wi = w(£)£ . It turns out that there is a closed form solution for the ^-integral, 

/>oo 

Wi = \ d£w(0Z l = r(Z + 1) ((I + 1) for photons, (B3) 
Jo 

where T is the Gamma function and £ is the Riemann zeta function. Some numerical values are given here: 

Wi » 1.645, W 2 ~ 2.404, W 3 w 6.494, W A » 24.89, W 5 « 122.1, W 6 ~ 726.0, etc, forphotons. (B4) 
Therefore, with our shorthands, 

Totally I indices 

'• / q \ 1+2 p 

R o ■ ■ ■ wife • • • i m = |£ ( _ \ dn [(eg Wi +1 a + C? Wi+ 2 a + C%W l+3 a 00 + • • •) n n n 12 ■ ■ ■ n 1 '" + 
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(c{Wi +2 a h + C*W l+3 a 0jl + C\Wi+ A a mi + ■ • >jn^n 12 ■ ■ • n^n^ + 
(c%W l+s a jlh + ClWi +A a 0jxh + ClWi +5 a 00jlj2 + ■ ■ -)n ll n i2 • • • n im n jl n j2 + 
(c$Wt +4l a jlj2js + ClWi +5 a 0jlhj3 + C^ +6 a oiu 2 i 3 + • ■ ^jn ll n 12 ■ ■ ■ n im n h n h n h + 

„ oo oo 

J n s - 



/ n \ '+2 c oo oo 

9sC f f 
h 3 



P ® 1 P Totally q indices 



\ 1+2 oo oo p 

9 -$[- r ) Y^1 ( V V <-;- "u Ohh ,. I '/<->»'''" ■■■„■„• n? 2 (B5) 



=0 q=p 

Totally q indices 



where C P denotes binomial coefficient C P = q\/p\(q — p)l. 

The integrands in equation iB5\ are products of different components of unit vectors. We can employ|Thorn3{l980)'s orthog 
onal condition to evaluate it: 

If 1 

— / dfln n n 12 ■■■n tn = 

47r / n + 1 



drin n n i2 ■■■n z " = — — 7 S% l2 '" tn , (B6) 



where we have defined Sl^ 12 "' ln be the Thorne delta function of rank-n. It carries interesting symmetric properties to help us 
simplify the evacuation of the integral. The function vanishes when n is odd. For even n, it i s defined to be the completely 
symmetrized product of Rronecker delta functions. For completeness, we copy the formulas from Thorne ( 198CJ): 

(5*1*2 -in = 5 {hh . . . £7»-lj„) = — V" gjljk 2 5 3k 3 j ki . . . (ji*„_lJ*„ ( B7 ) 

1 (n-l)\\ ^ 

v y k 2 k 4 ---k n 

where 

• ^2 is summed from 2 to n; 

• &3 is the smallest integer not equal to 1 or k 2 \ 

• &4 is summed over all integers from 2 to n, not equal to k 2 or £3; 

• feg is the smallest integer not equal to 1 or k 2 or £3 or ^4; 



We list all the non-zero Thorne deltas up to eight indices below: 

c -1 cii -1 ciiii 1 xiijj ^ riiiiii i riiiif? ^ xiijjkk 1 

(Vr = 1; <3 T = 1; d T = 1, \ J = -; d T = 1, d T = -, d T JJ = — ; 

1 Q 1 

^iiiiiiii _ ^ ^minjj _ j_ fiimjjjj _ ^mijjkk _ J_ (B8) 

T 35 35 

Note that the Einstein summation convention is not applied in the above equations. The symbols i, j, k denotes indices that are 
not equal to each other. With the help of the Thorne delta function, we obtain the expression for all moments 

Totally I indices 

, " , / n\ 1+2 oo oo (~ipTX/ 

Ra ^-.. ai ^ RQ ... 0ili2 ... im= ^9sC j^j ^ Y, ^TT Q .° • • ' ■ • • frff 8 "- 1 -**-* ■ (B9) 

P Q P Totally q indices 

h The number of zeros in the subscripts of a>o--.Qj 1 j 2 ---j p increase as we sum over q. 

The first f ew m omentum moments have important physical meanings. The first moment R a is the radiative four-flow. Based 
on equation I 



R ° 'IS 



o / n\ 3 
S7TC 



(BIO) 



f - J (W 2 a + W 3 a + W ia0 Q + .•■)+- {w iajlh + 3W 5 a 0jlh + 6W 6 a 00jlj2 + • • •) 5% h + ■ 
Because a aia2 ... ai has trace-free spatial part, the above equation reduces to 

R° = ^(w 2 a + W 3 a + W 4 aoo + ---), (Bll) 
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where we have defined the length scale A = he/ '9, which is proportional to the photon mean separation (i.e., n -1 / 3 , it is a purely 
radiation properties and has nothing to do with the material-dependent photon mean free path). Similarly, the spatial components 
of the radiative four-flow take the form 



R' 



8ttc 
3A 3 ' 



W 3ai + 2W 4 a 0l + 3W 5 



(B12) 



It becomes clear now, although a temperature is not well defined in non-equilibrium relativistic statistic mechanics, the notion of 
mean photon separation always exists. It is possible to use this physical interpretation to choose 9. The second moment R a ° is 
the radiation stress-energy tensor. Its components take the form 



R oa - 
R lJ - 



8irc9 , 

■ -(W 3 a + W 4 a + W 5 a 00 



2W 5 a 0i + 3W 6 a aoi 



A 3 c 
Wire 9 ( 

—^3 -\W 5 aij + ZWeaoij + 6W 7 a 0Qi j 



^R 00 6^. 



For R a ^, 



R 0Q0 = 



R 00i = 

R 0ij = 
R ijk = 



8ttc 
1? 
8nc 
3A 3 

16nc 
15A 3 
16ttc (9 



-^j (w 4 a + W 5 a + Weaoo + ■ 

W 5 cii + 2W e a Ql + 3W 7 a Qi + ■ 
(Wedij + 3W 7 a 0i j + OWaaooij 



+ 



:R 



000 x ij 



35A 3 



W 7 a ijk + 4W s a ijk + WW 9 a oijk 



\{R 0K 5^ 



R 00 jgki 



R 



(B13) 
(B14) 
(B15) 

(B16) 
(B17) 
(B18) 
(B19) 



Higher order moments can easily be obtained by using equation JB9t . Fixing the energy scale 9 at each point x a , the moments 
are simply linear transform of Grad's coefficients. Hence, we can solve for the coefficients by applying an inverse linear transfer. 
Although the zeroth moment R is never used in radiative transfer, for completeness, we give its expression here 



R = c J^lf = ^(w 1 a + W 2 a + W 3 a 00 + ---y 

C. GENERAL FORM OF THE MOMENT EXTINCTION TERMS 



(B20) 



The right hand side of the radiative transfer equation has the same mathematical form as the linearize d collision term of the 
Bhatnagar-Gross-Krook model (see standard textbook such as lHuandl!987l : ICercignani & Kremeril2002l) . We start deriving its 
moments by defining 



G a ia2 - ai = c f±P x r. _ s U«i p a 
J P 



■■■jar 



1+2 



d£dnx(f - s)£ l+1 n ai n a2 



(CI) 



where s denotes the source distribution function and x is the Lorentz invariant extension coefficient. Similar to the appendix iBl 
we define £ = hv/ksT and replace the integral over v by the integral over £. We assume thermal radiation so 



9s/h 3 



e 

S _ x e P a u a /k B T _ I 



H»<0, 



(C2) 



where we further replace U a by the fluid four-velocity u a and 9 by the fluid temperature k^,T for u>(£). Substituting iGradf s 
series for /, the extinction terms become 



G 



¥ irf) 1+2 J ^ dnxwi +! ( a - 1 + ^ nPl + fa^ 1 ^ + •••)' 



(C3) 



The zeroth and first moments of the collision term have clear physical meanings. When I = 0, G is the photon extinction 
rate; while for I — 1, G a is the radiative four-force. Higher order moment such as G af3 are simply called the extinction term of 
the third moment because of the balance equations. We evaluate (jf" 1 " 2 '""' in the fluid comoving frame so that u a has only the 
temporal component. 

The Lorentz invariant extension coefficient x = (hv/c)xv depends on the radiative process. Consider free-free, bound-free, 
or electron scatter, we notice that each of them is proportional to some power of temperature and frequency, \v ^ T^v^ ', in 
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standard approximations. Therefore, we assume for the general form, 

x ^im(hLY + *-\ (C4) 



A \m c c 2 _ 

where A = hc/k^T and r(£) is some dimensionless function in £. The momentum moment of the collision term due to a 
particular radiative process becomes 



G - A 4 



!t ' C (^F) (S)^ 1 J d ^ dnT (0w l+ 2(a-lH^y i +ea M2 n^ + ---)n^n^ ■ ■ ■ n*' . (C5) 

The expression is very similar to equation (IBU . We can further define 

b = a-l, bp 1 =ap 1 , & ft( g 2 = ap^, ... (C6) 

and the shorthand 

/>oo 

X l = \ dir{i)w{i)i l . (C7) 



Note that I is not necessary an integer in the above definition. Fortunately, the values of X\ are well define as soon as \v does not 
grow exponentially. Follow the same procedures in appendix|B] we can easily deduce the following expression 

Totally I indices 

A 4 V c J \m e c 2 ) ^ ^ p ™ + p + l ° " * 



Totally n indices 



The zeroth order moment is the photon extinction rate. It is simply 



G 



8nc ( k B T ^ 



A 4 \ m e c 2 

The radiative four-force takes the form 



)<p-t-y> — ± 
[X 2 b + X 3 b + X4&00 + •••)• (C9) 



3A 4 c V m c c 
For the second moment G a ^, we have 

&ttc fk B T\ 2 /fer^ -1 



G0= 8vrcfcBT/fcBTA - (j^ + ^ + j^ + ...), (C10 ) 
A" 1 c \m c c z J V / 

Gi = gcfc^T ^M_y +V ' 1 + 2X 5 b 0i + 3X 6 b 00t + •••). (Cll) 



G00 = ^c^j (X 4& + X 5 6o + X 6 6oo + -..), (C12) 

GOi = w- 1 ^ + 2XebQj + 3X7&ow + ^ (C13) 



(^Y {S)^ • 3X 7 & 0ij ■ 6X 8 b , ....)■ (C14) 



Other higher order moments can be obtained by using equation ( IC8b - If there are more than one radiative process in the problem, 
we can compute the extinction terms for each process and then sum over the results. 



D. SENSITIVITY OF FIDUCIAL VELOCITY 



We consider a very simple radiative transfer problem. Assuming there are only two streams of photons moving along ic 1 -axis 
in opposite directions, the photon distribution function is 



/ = || \n+5{v - v+)8 z {n l ~ 5 n ) + n-5{v - v-)5 3 {n l + S l1 ) 



(Dl) 



where n± and v± are the density and frequency for the two streams, 5{v) and S 3 (n l ) are the one- and three-dimensional Dirac 
delta functions. It is trivial to evaluate its momentum moments. The non-vanishing terms are 



I indices I indices 

? o~^"oo „ („ + p0---01 „ („ „ ,,1+1} 



i?°'-' 00 oc (n + ^ +1 +nV +1 ) , R° ' " ' 01 oc - nV +1 ) , . . . (D2) 
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For an Z-th order decomposition, we choose the fiducial reference frame so that the Z-th order flux i? " 01 vanish. This is equivalent 
to solve for (3 = U 1 /c in the equation 



;+i / , v 1+1 




n+ [ u+ ]Jh^ 1 =n -i^/^i ■ (D3) 

Suppose Pi is a solution to the above equation. To see how the fiducial frame depends on the distribution function, we suppose 
there is a small change to the photon number density n±. The corresponding change in the fiducial velocity is related to the 
derivative 

df3 | 

' ~"-2n±(l + 0' 



dn± 



= ± 1 Ji (D4) 



In the limit I — > oo, the derivative d(3/dn± — > 0, which recovers the linear closure. Although the general situation is more 
complicated, it is sensible to conjecture that the fiducial velocity obtained from higher order decompositions are less sensitive to 
the distribution function. Hence, the non-linear Lorentz transforms should introduce less unphysical photon self-interaction. 
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